arXiv: 1502.05449v4 [astro-ph.EP] 30 Apr 2015 


Draft version May 1, 2015 

Preprint typeset using DT^jX style emulateapj v. 5/2/11 


SPACING OF KEPLER PLANETS: SCULPTING BY DYNAMICAL INSTABILITY 

Bonan Pu (jftfilS) & Yanqin Wu (KSJA) 

Department of Astronomy & Astrophysics, 50 St. George Street, University of Toronto, Canada 

Draft version May 1, 2015 

ABSTRACT 

We study the orbital architecture of multi-planet systems detected by the Kepler transit mission 
using N-body simulations, focussing on the orbital spacing between adjacent planets in systems show¬ 
ing four or more transiting planets. We find that the observed spacings are tightly clustered around 
12 mutual Hill radii, when transit geometry and sensitivity limits are accounted for. In comparison, 
dynamical integrations reveal that the minimum spacing required for systems of similar masses to 
survive dynamical instability for as long as a billion years is, ~ 10 if all orbits are circular and copla- 
nar, and ~ 12 if planetary orbits have eccentricities ~ 0.02 (a value suggested by studies of planet 
transit-time-variations). This apparent coincidence, between the observed spacing and the theoretical 
stability threshold, leads us to propose that typical planetary systems were formed with even tighter 
spacing, but most, except for the widest ones, have undergone dynamical instability, and are pared 
down to a more anemic version of their former selves, with fewer planets and larger spacings. So while 
the high multiple systems (five or more transiting planets) are primordial systems that remain stable, 
the single or double planetary systems, abundantly discovered by the Kepler mission, may be the 
descendants of more closely packed high multiple systems. If this hypothesis is correct, we infer that 
the formation environment of Kepler systems should be more dissipative than that of the terrestrial 
planets. 


1. INTRODUCTION 


NASA’s Kepler mission has discovered that many of 


multiple planets (Kepler multis, 

Borucki et al. 2010 

Lis- 

sauer et al. 2011a Batalha et al. 

2013). Planets detected 

in these systems likely possess 

ow mutual inclinations 

(e.g. Lissauer et al. 2011b Tremaine & Dong 2012 

Fang 

& iVlargot 2012 

Fabrycky et al 

2014|), and are spaced 


planets (4-6), the spacing is so tight that there may be 
of order ~ 10 planets fitted inside the orbit of Mercury. 

Systems where a single planet transit (Kepler single) 
and where multiple planets transit (Kepler multis) may 
share the same intrinsic architecture but only differ in 
transit geometry. However, a few statistical studies have 
shown that this is not the case. There appears to be 
an excess in the number of singles over what one ex¬ 


pects by extrapolating from the Kepler multis (Lissauer 
et al. 2011a; Jo hansen et a l. 2012[ |Weissbein et al.pOKf 


Ballard & Johnson 2014|), suggesting that many of the 

singles are intrinsically singles. Even the multis may 
contain different sorts of systems. A direct probe of 
the underlying popula tion i s provided by transit-time- 
variations (TTV). |Xie et al. (2014) probed the presence 
of close-companions (transiting or not) around Kepler 
candidates using TTV signals. They found a striking 
correlation between the fraction of planets that show de¬ 
tectable TTVs and their transit multiplicity. Systems 
where four or more planets transit enjoy four times higher 
TTV fraction than Kepler singles do, and about twice as 
high as those with two or three transiting planets. This 
led them to propose that there are at least two differ¬ 
ent classes of Kepler systems, one closely packed and 
one sparsely populated (also see Weissbein et al. 2012). 
Besides from planet spacing, orbital eccentricity and in¬ 


clinations also differ between the two classes. Compared 
to Kepler multis, the transit durations of Kepler singles 
appear to be longer, requiring their orbits to have higher 
eccentricities (Wu, unpublished); study of stellar obliq¬ 
uity suggests that Kepler singles may lie on orbits more 


misaligned w ith the stellar spin than the multis (Mor¬ 
ton & Winn 1 2014). Both lend to the picture that Kepler 
multis are dynamically cold systems, while singles are 
dynamically hotter. 

What co uld lyave produced such a dichotomy? |Jo-| 


hansen et al. (2012) explored the possibility of dynamica 
instability. They boosted the masses of observed triple 
planetary systems by about a factor of 10 and found that 
these systems could become destabilized within 10 Bil¬ 
lion years. This led them to hypothesize that Kepler sin¬ 
gles may be the aftermath of dynamical instability. In 
this study, we provide much more definite arguments to 
this hypothesis by investigating the spacing distribution 
among the closely packed multi systems (4 planets and 
more). We find that, adopting realistic planetary masses 
and ages, Kepler high multiple systems are packed close 
to the boundary of stabilityJJ This proximity has to be 
explained. Unless the formation mechanism is somehow 
capable of producing this proximity, the most likely ex¬ 
planation is that the primordial systems had a range of 
spacing, but the ones that were packed too closely have 
been sculpted away by the subsequent billions of years 
of evolution. In this hypothesis, Kepler singles and lower 
multiples are the descendants of Kepler multis. If this 
true, we may be able to recover the primordial spacing 
of planetary systems. 


1 In contrast, the lower multiple systems are much more l oosely 
packe d (including the triple systems studied by | Johansen et al.| 
[ 2012 ). 


























































2 


2. HOW TIGHT CAN YOU PACK PLANETS? - 
THEORETICAL EXPECTATIONS 

2.1. Previous Works 

Here, we set the scene by reviewing some previous re¬ 
sults on stability in multiple planetary systems. 

Johannes Kepler (1571-1630), after whom the Kepler 
mission is named, initially believed that the spacing of 
planets in the Solar system is limited by packing perfect 
solids (polygons with identical faces). He discarded this 
idea after developing his famous Kepler laws. However, 
the question still remains, how tight can you pack planets 
together without affecting dynamical instability. Over 
the past few decades, some theoretical progress have been 
made in answering this question, and in some special 
cases the answers are known exactly. 

One usually measure spacing between every two plan¬ 
ets in a unit called the mutual Hill sphere 


Rh = 


i / i \V 3 

cii + a 2 / mi + TO2 \ 


V 3 M* 


( 1 ) 


where dj, rrii are the semi-major axis and mass of the 
two planets, and M* the mass of the central star. In this 
unit, planet spacing can be expressed as a dimensionless 
number I\ 

K =^ ( 2 ) 

Kh 

We first focus only on orbits that are circ ular and 
coplanar. A well known result, derived by |Wisdoni| 
(1980) based on overlapping of first order mean-motion 
resonances, is that a test particle (m 2 = 0) perturbed by 
a planet (mi, aq) is subject to orbital chaos if it lies too 
close to the planet, 


|Q2 - ai| 

(02 + ai )/2 


< 1.49/j, 2 / 7 . 


( 3 ) 


where /i = mi + m - is the planet mass ratio and the nu¬ 


merical coefficient is adopted from Duncan et al. (1989). 
In K value, the critical spacing is 


A'crit = 3.46 


P 


4.5 x 10 


—5 


-(A) 


( 4 ) 


This result is little modified i f the test particle is in¬ 
stead substituted by a planet (Deck et al. 2013). The 
dependence of AT cr i t on p (albeit weak) indicates that 
Hill sphere is not necessarily the best unit to measure 
stability. However, we adopt it here for convenience. 
There is also an important result on Hill stability (stabil¬ 
ity agains t orbit crossing), popularized in astronomy by 
Gladman (1993). As a result of total energy and angular 


momentum conservation, two initially circular, coplanar 
planets can avoid close encounters at all times if they are 
separated by more than A' cr i t = 3.46. These two criteria 
meet at p = 4.5 x 1 0~ 5 and their physical connection is 
discussed in Deck et al. (2013). 

For systems where there are more than 2 planets, there 
is yet no theoretical understanding for the o rigin of dy^ 
namical instability (see attempts by, e .g. Chambers et al. 
1996 Zhou et al. 2007 Quillen 20111. Whatever infor¬ 
mation we have on how tight one can pack planets in 
these cases are empirical results from numerical simula- 
tions (Chambers et al.(l996; Yoshinaga et al.|1999; Zhou 


et al.|2007|lChatterjee et al.|2008[ [Smith fc Lissauerj2009j 

Funk et ahj2010). These we summarize below. 


Generally speaking, systems with more planets (N) re¬ 
quire larger K spacing to stay stable. To compare the 
various numerical studies, we scale the physical timescale 


tern, r = t/Ti (also see 

Funk et al. 

2010 

). We also sup- 

press error-bars in the fitting resu 

Its but caution that 


there are intrinsically large scatters. 

• dependence on r: Two-planet systems can remain 
indefinitely stable as long as they are above a criti¬ 
cal spacing. This may not hold for large-N systems. 
It appears that a larger spacing does not guarantee 
eternal peace but only delay s the time of instability 
(we return to this point in f 2.31. Numerically, the 


critical spacing at which an average system under¬ 
go e s their first encounter at time r can be modelled 


as (Chambers et al. 1996 
Funk et al.||2010 p 


Smith & Lissauer 2009 


A„-it = a + b log r , 


( 5 ) 


though the alternative scalings log K = a + b log r 
and K = c l og[b(l ogr — q)1 have been adopted by 
Zhou et al.| (l2007|) and |Chatterjee et al.| (|2008|) re¬ 
spectively. The numerical coefficients a and b de¬ 
pend on system parameters like the number of plan¬ 
ets, their masses, and orbital eccentricities and in¬ 
clinations. Two results most relevant to our case 
are 


AT crit « (1.33 ± 0.3) logr + (0.0 ± 2.0). 


( 6 ) 


for systems of 8 planets all with p = 5x 10~ 5 (Funk| 
et ai'| |2010), and 


AT cr it ~ k>gT + 1.7 


( 7 ) 


for 5 planet systems with p = 3 x 10 -6 (Smith 
fe Lissauer||2009 Zhou et ah1|2007 ). Both scalings 
are obtained for circular, coplanar systems with a 
uniform planetary mass and a uniform K-spacing. 

dependence on N: The most com prehensive stud y 
to date on this issue is that by Funk et al. (2010). 
They found, for r = 10 8 ' 5 , and p = 5 x UP ', the 
critical separation for a system of circular, coplanar 
planets goes as 


AT 


crit 


4 + (N — 2) = 2 + N if A? < 8 
10 if N > 8. 


( 8 ) 


This results roughly recovers equation (HI) when 
N = 2. The initial rise of AT cr i t with N )is intu¬ 
itively understood as the presence of more plane¬ 
tary perturbers exerts more strain on the delicate 
dynamical systems. The subsequent flattening of 
AT cr it with N is, however, surprising and not un¬ 
derstood. One guess is that dynamical instability 
is most affected by close neighbours, and the pres¬ 
ence of far-away neighbours, once beyond a certain 
separation, do not matter. For instance, a system 
of 8 Neptunes will space a period range of « 10. 
With such a large spacing, newly added planets 
have no chance of participating in any MMRs with 
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the pre-existing planets, and therefore should not 
adversely affect the stability of the system. 


• dependence on fi: It is found that less massive plan¬ 
ets, with their smaller Hill spheres, need b e sepa¬ 
rated by slightly hig her values of K. B oth Cham- 
bers et al.| dl996 ) and|Zhou et al.|(|2007|) found that, 
for every decade of drop in p, (within the range of 
p = 10~ 9 to p = 10” 3 ), the critical separation K 
should increase by 0.5 — 1. 

• Eccentricity/Inclination: If the planetary orbits 
are dynamically hot, mor e elbow room is neede d 
to maintain stabil ity (e.g. Yoshinaga et al. 1999). 


Zhou et al. (2007) showed that, at low eccentric¬ 
ities, increasing the average eccentricity by 0.01 
raises K cr i t by ~ 1. We know of no study that ex¬ 
plores the impact of mutual inclination separately, 
so we provide some illustrations of this below. 


Returning to our problem at hand - a system of 8 flat, 
circular planets, with masses of ~ 5M®, and an inner 
orbit of 0.1 AU, may remain stable for over 1 Gyrs (or 
r > io 10 - 5 ) , if K ~ 14 ± 5 (extrapolating from equation 
[6]) . This valueis broadly consistent with the spacing ob¬ 
served in tight Kepler systems. However, a more refined 
approach is required. 


2.2. Our numerical explorations 

To compare theoretical results to Kepler systems, we 
need to overcome a few restrictions of previous studies. 
First, in almost all previous numerical simulations, plan¬ 
ets in a given systems are assigned identical masses and 
identical spacing, in contrast to realistic systems. Inho¬ 
mogeneities in mass and spacing may impact results sig¬ 
nificantly, e.g., they will introduce large scatter in stabil¬ 
ity lifetimes, among systems that have comparable bulk 
parameters. It is therefore necessary to study system sta¬ 
bility in a statistical sense - previous investigations have 
studied only a small number of cases (typically fewer 
than a hundred systems). We accomplish this by inte¬ 
grating a much larger number of systems (~ 10 4 ). In 
addition, we wish to study the impacts of small eccen¬ 
tricity and inclination. Lastly, our simulations run for a 
factor of 10 longer than previous works, up to 10 8 years. 
Although still shorter than the typical system ages (a 
Giga-year), they provide a better indicator for long-term 
stability. 


as we know the stabilit y thr eshold only depends weakly 
on the planet masses ((2.1). The inner most planet is 
placed at 0.1 AU. We then determine the spacing, in 
unit of mutual Hill sphere, for an adjacent pair of plan¬ 
ets by drawing from a normal distribution with mean 


K mean and variance <ta', 


| _ (K — Kmean) 

P(K)dK = - =e 

(JkV 27T 


dK. 


(9) 


The initial arguments of periapsis and the mean motions 
are drawn from a uniform distribution in [0, 27 t] . 


We adopt the Wisdom-Holman integration scheme 

Wisdom & Holman 


1991 

), as implemented in the 
. E. Kaufmann, an improved 

SWIFTER package by 


Javic 


version of the original SWIFT package by Levison and 
Duncan. The post-Newtonian correction is added as an 
additional potential of the form Uqr = —3 (GM+/cr) 2 , 
and we adopt a time-step of 10 -3 years. Integrations 
are discontinued either when two planets approached one 
another within one mutual Hill radius (and therefore are 
destined for a close encounter), or when a preset time 
limit (usually 10 7 years) is reached. We record the sta¬ 
bility time of each system, and analyze how it depends on 
various parameters. Of particular interest is the survival 
probability at a given time. 


2.2.2. Results for Flat, Circular Systems 

We integrate, up to a time of 10' years, a total of 
8,000 circular, coplanar systems (cr e = <7,; = 0) with 
6.0 < K m ean < 13.0 and 0 < cjk < 4.0. In addi¬ 
tion, we simulate 120 such systems, within the range 
6.0 < K mean < 13.0 and 3.25 < cr K < 3.75, up to a 
time of 10 8 years. 

We report the survival probability, the fraction of sys¬ 
tems that remain stable up to a time t, as a function 
of planet spacing. Overall, a wider K me an leads to a 
higher survival probability (Fig. [I), however, the dis¬ 
persion in K (<jk) also impacts stability. For the same 
Kmean , systems with a greater K dispersion amongst 
its pairs are less stable. In other words, a pair of planets 
with a smaller-than-average I\ have a destabilizing effect 
on the whole system that out-weighs any corresponding 
stabilizing effect afforded by a planet pair with a larger- 
than-average K. 

We find, empirically, that our results can be succinctly 
summarized using an effective spacing, 

K e s = Kmean 0.5fX/c ■ (10) 


2.2.1. Numerical Setup 

We generate artificial planetary systems as follows: 
each system contains 7 planets orbiting around a 1 solar 
mass star, on coplanar and circular orbits. The number 
7 is chosen for a number of reasons: this is the largest 
num ber of tran siting planets in a multiple system (e.g. 
Lissauer et ah_ 2014) ; it is close enough to the plateau 
in equation ® , the results may be considered valid for 
all systems with more than 7 planets; it is computation¬ 
ally feasible. The mass of the host star is taken to be 
M* = M 0 , while the mass of the planets are drawn from 
a uniform distribution between 3.0 Me > m > 9.0 Me, 


We arrive at this combination by noticing that, in our 
simulations, increasing cjk by 1.0 has the same effect 
on stability as decreasing K mean by 0.5. So for systems 
with a varied spacing, K e g acts as an equivalence of K in 
systems with a uniform sp acing. Such a combinat ion can 
also explain the results in Chambers et al. (1996) where 
they made some preliminary explorations on the impact 
of a Jv-dispersion. 

We can now fit the survival probability reported in 
Fig. [ 2 ] as a function of K e g (the functional form being 
somewhat arbitrary) 


s(K e g, t ) = 1 - g- [°.5iC e ff-0.6—0.35 log t ] 2 


( 11 ) 


systems (Wu & Lithwick 

2013 

Hadden & Lithwick|2014 for K e g > 1.2 + 0.7 log r and 0 otherwise. Here the di- 

Weiss & Marcy 

20MfT r J 

L'he latter assumption is robust mensionless time r = t/1 i, and 1\ is the orbital period 
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of the inner most planet (~ 12 days in our case). We 
now define a threshold spacing, K 50 (t), as the smallest 
K e ff necessary to achieve a survival probability of 50% at 
time t. This spacing has a simi l ar fun c tional form to that 
presen t ed in|Chambers et al.| (fl996); Smith & Lissauer 
(2009); Funk et ai. (|2010|), but with different numerical 
coefficients, 


K 50 (t) » 0.7 log t + 2.87. 


( 12 ) 


We compare our critical spacing with those from previ¬ 
ous works in Fig. [3] Our results, with a spread in K and 
in planetary masses, have to be characterized by A' 50 . 
This differs from the definition of /\ cr it for other groups, 
which typically adopt a uniform I\ and mass. Ho wever, 
results br oadly agree within error-bars (as given by Funk 
et al. 2010) but ours tend to fall below. In other words, at 
the same spacing (measured by K e g in our cases), inho¬ 
mogeneous systems (with mass and spacing dispersions) 
remain stable slightly longer than homogeneous systems 
do. We note that homogenous systems, with their uni¬ 
form period ratios, tend to be strongly a ffected by the 
presence of mean-motio n resonances (Smith & Lissauer 
2009 Funk et al. 2010 see, e.g.). This may partially 
explain the difference. A more important difference be¬ 
tween homogenous and inhomogeneous systems is that 
the former may reac h eternal stability at a f inite spac¬ 
ing (as suggested by|Smith fe Lissauer 200 9), while we 
do not observe this for the latter (see ' ]2W|) From now 
on, we consider that our results (equation|l2[) supercedes 

equations B :0 - . 

Our scaling is valid for 4 orders of magnitude in r (r = 
10 5 - 5 to 10 9 - 5 , Fig. (2J). In subsequent sessions, we assume 
that it remains valicl for another one order of magnitude, 
to r = 10 10 5 (or t = 1 Gyrs). At this epoch, the critical 
spacing is K 50 ~ 10.2. 


2.2.3. Results for Eccentric, Inclined Systems 

Here, we explore quantitatively how eccentricity and 
inclination affect the critical spacing. We generate 12,000 
artificial 7-planet systems with a normal distribution in 
spacing, 6.0 < K mean < 20.0, 3.0 < ctk < 4.0, and 
a Rayleigh distribution in eccentricity and inclination. 
The scale factors are a e and <Ti„ c , e.g., 

e AT 

P(e)de=^e 2 ^de. (13) 

We explore the range a e £ [0,0.05], and cr inc £ [0,0.16] 
(in radian, corresponding to a range of [0°,9°]). The 
initial longitudes of ascending nodes are randomly dis¬ 
tributed between 0 and 27r. 

Given the extra dimension in orbital parameters, it is 
only feasible to integrate these systems up to a time of 
one million years (r = 10 7 °). The survival probability 
for different types of orbits are shown in Fig. [4] 

Individually, mutual inclination and initial eccentricity 
each has a strong destabilizing effect on stability. For the 
range of parameters that we explore, the critical spacing 
can be characterized as 

K 50 (a e ,a inc ) « K 5 q(0,0) + (^_) + (^) , (14) 


where X 5 o( 0 , 0 ) is the corresponding value for circular, 
coplanar systems for the same r (eq. 12). An eccentricity 


Planets per SMA decade 



Fig. 1.— Survival probability, measured at time t (with t rang¬ 
ing from 10 4 to 10 8 yrs), for systems with different K mean . Each 
point represents results from some 500 systems, with its horizontal 
error bars representing the width in Kmean and vertical error bars 
accounting for Poisson error in numbers. The solid lines show the 
aggregated results for a distribution of cr^. The latter is taken to 
be a Gaussian distribution with a mean of 2 and a dispersion of 
~ 1, with the exception of the black curve where 3.25 < ctk < 3.75. 
It appears that K mean is a reasonable indicator for system stabil¬ 
ity. However, system stability also depends on the dispersion in 
K. As the two curves (dashed and dotted, both for t = 10 7 yrs 
but each with a different ctk) show, systems with larger dispersion 
tend to be more unstable. So simply measuring the mean does not 
reveal the whole picture. 



describe the nu meric al results succinctly. Different solid lines rep¬ 
resent equation ( |llj ), evaluated at their respective times. Dashed 
curve shows the extrapolation of equation ( |ll| ) to t = 10 9 yrs (or 
T=10 10 - 5 ). 


dispersion of a mere one percent requires the spacing to 
expand by one Hill sphere, for /i ~ 10 -5 . In relative 
terms, mutual inclination is less destructive to stability, 
by about a factor of 4. 

One can also scale the above result in a natural unit 
of the problem, the Hill eccentricity, en = Rh/u- For 
our parameters, [fii ~ ~ 2 x 10 -5 ), en ~ 0.024. The 
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Fig. 3.— The critical K spacing required to maintain stability 
over time r (x-ax i s, i n logarithm), as obtained by different groups. 
The |Funk et al.| |2010 ) results (in red, equation |6| with vertical 
error bars) are for ll = 5x 10 5 ; the |Smith &; Lissauer | < |2009| ) results 
(in blue, equationm) applies to u = 3x 10"°; we insert /x = 2x 10 -5 
into the fitting formula of|Zhou et al. (2007|), and modify their 
expression to suit present notation, to obtain the magenta line. 
All these simulations are of uniform mass and uniform spacing and 
therefore are not statistical in nature. Our results, K $o (eq. |12| ), 
are plotted as a black curve. The horizontal extent of each curve 
also delineates the respective time range probed by different studies 
- the present study probes up tor = 10 8 ' 5 , ten times longer than 
previous ones. Our critical spacing are broadly consistent, though 
somewhat smaller, than those of previous studies. 

expression for critical spacing, re-cast in these units, is, 


K 50 (a e , dine) ~ K 50 (0, 0) + 0.4 x 


eH 


4.0en 


(15) 

This leads to a surprisingly geometrical interpretation 
for our result. For every unit of eccentricity rise (in unit 
of en), the average spacing between two adjacent plan¬ 
ets (K ) have to rise by 2.5 Rh in order to maintain sta¬ 
bility; however, such a rise in eccentricity also reduces 
the closest approach (measured from the apo-apsis of 
the inner planet to the periapsis of the outer one) by 
~ 2ena = 2i?n- So a system of mildly eccentric planets 
behaves similarly as a circular system that have the same 
closest approach. We do not have a ready explanation 
for this result. 

In the following, we assume that the survival stability 
obtained at these shorter integrations can be extrapo¬ 
lated to longer times, in the form of Equation (12). 


2.3. Discussion: Extrapolating to 1 Gyrs? 

We do not yet understand the origin for the scaling in 
equation (12). So it may be prudent to ask: can the scal¬ 


ing be extrapolated to as long as one billion years? For 
instance, in two-planet systems, it is known that beyond 
certain spacing the system can remain eternally stable. 
If high-N planetary system exhibits the same behaviour, 
equation ( |12[ ) will break down at a certain spacing. 

Lacking analytical under standings, we turn to num er- 
ical simulations for insight. Smith & Lissauer (2009) ex¬ 
perimented with equal-mass, constant-spacing 3-planet 
(or 5-planet) systems and found that as planet spacing 
reaches K ~ 7 (~ 8 for 5-planet), the stability time, 
clocked at 10 7 dynamical time at this point, starts to 
lift off and increases with K much more steeply than Eq. 



Fig. 4. — The critical spacing a; 
entricity, at time 1 Myrs (r = 10 7 ' 5 ). 


as a function of orbital ec¬ 
centricity, at time 1 Myrs (r = 10 /0 ). Different curves stand for 
different mu tual inclinat i ons. F or comparison, we overplot in dots 
the result of Zhou et al.| ( |2007j ) for coplanar systems. 


(12). It appears that this equation is violated beyond 10' 
dynamical times, and indefinite stability at large spacing 
appears possible. 

However, our results here suggest the opposite. We 
have integrated non-equal-mass, non-uniform spacing 7- 
planet systems for up to 10 9 5 dynamical times, while still 
observing that Equation (12) holds well. There is no un¬ 
usual transition that occurs at 10 7 dynamical times. We 
suspect that our heterogeneous set-up (a dispersion of 
mass and spacing) opens up a more complex phase-space. 
Such dynamical system may not exhibit a single charac¬ 
teristic timescale, as the more restrictive cases in Smith 


& Lissauer (2009). In which case, it may be reasonable 
to extrapolate Eq. © by one order of magnitude, to 1 
Gyrs (10 10 5 dynamical time). 

2.4. Discussion: Effects of Mean Motion Resonance 


Previous studies (Chambers et al. 1996 Funk et al. 


20101 have identified that mean motion resonances tend 
to be disruptive for dynamical stability. This i s also 
demonstrated in more detail in |Mahajan & Wu| (2014) 
where they showed that, when the positions tor Kepfer- 
11 planets are shifted away from their current values, 
they may encounter mean motion resonances and become 
destabilized. Here, using the 8000 circular, coplanar sys¬ 
tems that we have produced, we perform a statistical in¬ 
vestigation on the effect of resonances on stability (Fig. 

li¬ 
lt appears that, by 10 7 yrs, systems that contain near 
resonant pairs are preferentially destroyed, compared to 
neighbouring systems. This is an order unity effect for 
first-order resonances, and weaker, albeit still visible, for 
second order resonances. Similar carving is also seen in 
our eccentric, inclined ensembles. 

Amongst Kepler’s multi-planet systems, there has been 
an observed tendency for pla nets to ‘pile u p (just wide of 
mean motion resonances (Lissauer et al.|2011a). A num¬ 
ber of mechanisms have been pro posed for this asymme¬ 
try, including resonant repulsion (Lit hwick & Wu |2012 
Batygin & Morbidelli 2013 Delisfe et al. 2012| ), stochastic 
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1.2 1.4 1.6 1.8 2 

Period Ratio 


Fig. 5.— The upper panel depicts the survival rate vs. period ra¬ 
tios for all pairs of planets in circular, coplanar simulations. Here, 
the survival rate is the ratio of stable pairs (at 10 7 yrs) The lo¬ 
cations of the first order (dashed red lines) and the second order 
resonances (dashed blue lines) coincide with reduced survival rates, 
showing that mean-mot ion resonances tend to de-st abilize systems. 
The little inset at the top-right focuses on regions of width d=2.5% 
around 3 first order MMRs (2:1, 3:2, 4:3). There is a definite asym¬ 
metry: systems on the near side of resonances are preferentially 
depleted relative to those on the far side. The lower panel is the 
observed period ratios among all Kepler systems. These also show 
a relative deficit of systems near these resonances in the near-side. 
The similarit y betw een the two is intriguing. 


migration (Rein 2012), adiabatic mass accretion (Petro¬ 


vich et al.||201 3), and planetesimal interactions (|Ciiat- 
terjee fe Fordj|2015). We observe that (inset in Fig. [5]) 


tor circular, coplanar systems, systems with period ratios 
just inward of MMRs are more severely carved relative 
to those on the far side of resonances. The asymmetry 
is more pronounced in slightly eccentric, inclined sys¬ 
tems. This asymmetric carving may present another ex¬ 
planation for the observed systems. More investigation 
is planned in the future. 


3. ORBITAL SPACING OF KEPLER’S MULTI-PLANET 
SYSTEMS 

3.1. Sample Definition 

Observational data from the Kepler Space Telescope 
was retrieved from the NASA Exoplanet Archive in June 
19, 2014. All planet candidates with a KOI disposition of 
’’False Positive” were removed. We are only concerned 
with high-multiple systems ( N p > 4 planets) and our 
sample include 4 six-planet systems, 18 five-planet sys¬ 
tems and 55 four-planet systems. There are a total of 
257 adjacent pairs, with 92 pairs in 5 P and 6 P systems 
alone. 


3.2. The intrinsic K distribution 

To compare against theoretical results in fj2] we require 
the intrinsic A'-distribution (as opposed to the apparent 
one) of Kepler’s multi-systems. However, such a distribu¬ 
tion is not straightforward to obtain, due to three factors. 
First, the K parameter depends on planetary masses, 
which are not known for vast majority of planets. Sec¬ 
ondly, the presence of non-transiting (‘missing’) planets, 


Planets per SMA decade 

15 10 7 5 4 3 



Fig. 6. — The observed redistribution of high-multiple Kepler 
systems are shown as solid histograms for the 4P and 5P+6P 
systems. The mass-radius relation is assumed to be M rj 
3Mr ir (A'/Tt© ). According to KS test, the two histograms are sta¬ 
tistically different, with the 5P+ systems more tightly spaced. 
The dotted curve shows the 5P+ systems when we adopt instead 
M = Mca( R/Ret .) 2 ' 0 ®. a scaling obtained for Solar system planets 
jl.issaner et al. 2011b). This leads to a more diffuse K-distribution. 


either because they are too small to be detectable, or 
because they are unfortunately inclined relative to our 
line-of-sight, will cause the observed K to be artificially 
inflated. Thirdly, one must account for selection effects. 
The KOI sample is severely in complete for planets with 
period outward of ~ 200 days (Petigura et al.|[20l3 Sil- 


burt et al. 2014). Systems that have tighter-than-average 
spacing will be preferentially picked up as high multi¬ 
ple systems, while those with larger spacings will not. 
This tends to skew the apparent AT-distribution toward 
smaller K values in these systems. 

We overcome the first hurdle by adopting the mass- 
radius relationship as obtained for planets that have mea¬ 
sured ma sses. Recent studies using both transit timing 
variation (Lithwick e t al.|201 2; Hadden fe Lit hwick 20141 
and radial velocity (| Weiss & Marcy 20I4|) techniques 
have shown that planet masses scales approximately lin¬ 
early with planet radii, or M/M© « 3.0A/A®. This ap¬ 
plies to planets with sizes of a few Earth radii, the vast 
majority of our sample, and recent studies show that this 
relation may even ext end to planets as small as 1.4A© 
(Dressing et al. 2014). This relation differs from the 
one tha t applies to bolar system planets |Lissauer et al.| 
( 2011b|k_ We plot the resulting apparent AT-distribution 
in Fig. [6] There appears to be a sharp rise in the number 
of pairs above K = 10, and a sharp drop beyond K ~ 22, 
with some 75% of pairs spaced between 10 > K > 22. 
In addition, the apparent A'-distribution for 4P systems 
and that of 5 P, 6 P systems are statistically different (KS 
test returns a value of 0.003), with the former peaking 
at K ~ 17 while the latter around K ~ 13. It appears 
that even 4P systems may be heterogeneous, containing 
both tight spacing systems and large spacing systems. 
This prompts us to focus our attention exclusively to 5 P 
and 6 P systems in subsequent sections. Lastly, we notice 
that our adopted mass-radius relation leads to a sharper 
peak in A'-distribution, when compared to other power 
laws. 
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To overcome the second and the third hurdles, it is 
necessary to “observe” synthetic populations to extract 
their apparent AT-distributions. We now describe our 
procedure in constructing and ’observing’ these synthetic 
populations. 

Each planetary system orbit around a central star that 
is Sun-like, and have planets whose radii are distributed 


based o n the deter mination in Petigura et al. (2013); Sil- 


burt et al. (2014), namely, with equal probability per 


logarithmic bin between 1 and 2.5i?®, and a linearly de¬ 
creasing probability between 2.5 to 15f?0. The mass of 
each planet is then determined using our default mass- 
radius expression. The orbital inclination, measured rel¬ 
ative to a fixed plane, satisfies a Rayleigh distribution 
with scale parameter <7i nc (eq. 13). Planets are placed 
in circular orbits, and all Euler angles are uniformly dis¬ 
tributed. The period of the inner-most planet in each 
system is drawn from a Rayleigh distribution with a scale 
parameter of 5 days. This distribution is a good repre¬ 
sentation for our sample of planetary systems, but our 
conclusion is not sensitive to this particular choice. To 
assign periods to the remaining planets, we assign a K 
value that satisfies Equation (IM with a mean K mean and 
a variance ax ■ Each system is assumed to extend out to 
well beyond the observing boundary. 

With these choices, our average system, with Pj n ~ 5 
days and m ~ 5 Me, contain roughly 7 planets (ranging 
from 6-10 planets) inward of 200 days. For the entire 
ensemble, the period distribution is flat in logarithmic 


space, consistent with observational determinations (JSil- 
burt et ah|2014 ). 


We then identify planets that transit (impact parame¬ 
ter smaller than the stellar radius). Not every planet that 
transits may be identified by the Kepler pipeline. We cor¬ 
rect for this by adopting the detection completeness as a 
function of pla net size a nd orbital period, as obtained in 
Fig. 2 of Silburt et al. (2014). This effectively limit the 
detection sphere to periods inward of ~ 200 days, and 
miss some small planets. 

Our procedure therefore accounts for both the missing 
planets and the selection bias. The A'-distribution from 
these ’detections’ are then compared against that in Fig. 

a We find that the best fit, using only 5P+ systems, is 
tained when A' mea n = 14, <Jk = 3.4. This corresponds 
to an effective spacing of K e g = 12.3. 

As is shown in Fig. [7J the ’observed’ A'-distribution is 
always broader than its intrinsic counterpart. The broad¬ 
ening is caused by intervening planets missed by transit 
surveys: our simulations show that for low inclinations 
{(Jin C ~ 2°), ~ 35% of the apparent 5P and 6P systems 
contain one or more missing planets hidden in-between 
the transiting planets. 

Moreover, the apparent tight spacing in high multiples 
is not a result of selection bias, namely, only systems that 
have intrinsically the tightest distribution can fit inside 
the detection sphere (~ 200 days) and be observed as a 
high-multiple system (see §41). 

We can translate the ATaistribution to a more intu¬ 
itive unit, the number of planets within a decade of 
semi-major axis, typically located from 0.1 to 1.0 AU 
for our case. This number is marked as ’planets per 
SMA decade’ in multiple figures. The Kepler 5P+ sys¬ 
tems contain between 6 to 10 planets per SMA decade, 
with a mean of ~ 8. Contrast this with only 4 planets 


Planets per SMA decade 



Fig. 7. — The underlying spacing distribution. The observed 
spacing, for 4P and 5P+ systems respectively, are shown as solid 
histograms, the green curve is our best-fit solution for the under¬ 
lying distribution of the 5P+ systems, and the blue curve its pro¬ 
jection into observable space by synthetic observations. Our best 
fit, Kmean = 14, (Tk = 3.4, for the 5P+ systems fail to account for 
some of the more widely spaced 4P pairs. Here, we adopt cr inc ~ 2° 
- the A'-distribution does not vary much when the inclination dis¬ 
persion is varied. 

per decade in the terrestrial region of the Solar system. 

Our determination of the intrinsic A'-distribution is not 
affected by our choice of inclination dispersion, perhaps 
a bit paradoxically. A higher dispersion means more 
planets in a given system escape detection in a transit 
survey, and this should have increased the average K 
value. However, this also disqualifies such a system from 
our ’high multiple’ sample. As a result, the intrinsic A'- 
distribution we obtain, by focussing on high multiples 
only, is not much affected. 


4. DISCUSSION 

Gyro-chronology studies, the dating of stars using their 
stellar rotation periods, reveal that stars in the Kepler 
field, at least ones that have measurable r otation peri¬ 


ods, have typical ages of a couple Gyrs (Walkowicz & 


Basri 2013 McQuillan et al. 2014). This motivates us 
to assume a default age of 1 Gyrs for all Kepler systems. 
This may be an underestimate in many cases. Older stars 
are much less magnetically active and therefore less vari¬ 
able. They also rotate with longer periods. As a result, 
they may not be selected for gyro-chronology studies and 
results from the latter may be biased towards younger 
stars. Fortunately, an older age simply strengthens the 
argument we present below. 

If all K eple r high multi systems are circular, coplanar, 
equation (12) predicts that 50% of the planetary systems 
can remain stable to beyond t = 10° yrs (r ~ 10 10 ’ 5 ) if 
A' e ff > 10.2. This is surprisingly close to our best fit for 
Kepler high multis, K e g = 12.3. 

The proximity is even more striking if one allows for 
just a small amount of eccentricity and non-coplanarity 
in th e sy s tems. Ado pting the measu rements by |Wu fc 
Lithwick (2013); Hadden fe Lithwick (20141 for TTV 
0.02, equatiofn[T4|) requires that AT 0 ff > 12.2. 


pairs, <j e 

This observation, that Kepler high multis are perched 
at the cliff of instability (see also Lovis et aL||2011 Fang 


& M argotp01 3), as we suggest below, is most likely the 
















































result of continuous sculpting. This suggestion is the 
main thesis of this work. 

Our intrinsic distribution exhibits a sharp fall-off at 
larger K (see Fig. f7|, similar to the observed population. 
We believe this faltoff is genuine, not a result of observa¬ 
tional bias. When synthetically observing populations of 
planets with a range of K me an and a k, we find that even 
systems with spacing as large as K mean = 17.5 (fixing 
( 7 A' = 3.5) can be observed as 5P+ only 40% less fre¬ 
quently than our best fit (A' mea n = 14). Their presence 
should therefore be detectable. The lack of such wide 
spacing system in the observed high multis is therefore 
likely genuine. 


4.1. Dynamical Scuplting 

The formation of Kepler planets themselves is not yet 


dTerquem & Papaloizou 2007; Ida & Lin 

2010 

Hansen 

& Murray 20126 IChiang &; Laughlin 

2013 

Gha 

tterjee & 

Tanj2014). However, it looks likely t 

lat t 

iey were made 


ence of hydrogen envel opes on these low-m ass planets, 
for bot h Neptunes (e.g. Borucki ct al. 2 0Tot and su per- 
201 .'{ ' Weiss & Marcy| 20 111. 


Earths (TWu & Lithwick 


How does such a formation time jibe with the orbital 
spacing of the observed systems? The nascent disk can 
weed out (and possibly reform) systems that are unstable 
within a few million years (/v e ff < 8.1, eq. 12), but it 
should not discriminate among systems that have spacing 
wider than this. On the other hand, while it is natural 
that only systems with spacing wide enough to be stable 
over their lifetimes can be observed, it is striking how the 
observed systems skirt the edge of stability. 

This chain of thought then leads us to propose the fol¬ 
lowing scenario. Planetary systems were initially formed 
with a range of K e g >8.1. This distribution is contin¬ 
uously sculpted by dynamical instability. So by 1 Gyrs, 
the remaining population will exhibit a sharp cut-off at 
the border of stability. We call this ’the sculpting hy¬ 
pothesis’. 

Under this hypothesis, it is possible to reconstruct the 
initial planet population, after we posit a number of as¬ 
sumptions, some more sound than others. We assume 
that the all systems have the same small AT-dispersion 
as we observe in today’s high multis, uk = 3.4. Since 
orbital spacing in a system likely reflects its formation en¬ 
vironment, we expect the spacings to be similar within 
a given system. We further assume that all systems at 
birth have the same dispersion in eccentricity as the cur¬ 
rently observed TTV sample, <r e ~ 0.02, and a compa¬ 
rable dispersion in inclination. To ’reverse’ the effect of 
carving, we augment the current planet population by 
three populations that are similar in number but shifted 
to a smaller K e g by 0.7 x log(f/lGyrs), for the three log¬ 
arithmic time intervals between t = 10 6 and t = fO 9 yrs. 
The value 0.7, according to our numerical experiments, 
is the range of spacing that is carved away by the passing 
of each logarithmic interval. And the initial t = 10 6 yrs 
is the putative lifetime of our nascent disks. 

The resulting initial population is illustrated in Fig. 
[sj In this model, by I Gyrs of age, roughly 75% of the 
original systems have been disrupted. This fraction could 
rise or drop depending on what we assume for the initial 


Planets per SMA decade 



Fig. 8. — The sculpting hypothesis. The solid histogram and 
the solid curves are those displayed in the bottom panel of Fig. [T] 
namely, the observed (histogram) and the underlying (green curve) 
spacing distributions for Kepler 5p and 6p systems. The dashed 
curves, from bottom to top, represent our reconstruction of the ear¬ 
lier spacing distributions at, 10 8 , 10 7 and 10 6 years, respectively. 
See text for details. The primordial population (posited to be the 
one at 10 6 yrs) are in average more tightly spaced, and contain 
roughly 4 times more systems than those remain today. 

A" e ff distribution. 

What happen to systems that experience dynamical 
instability, in the interval since their birth and today? 
Once two planets are excited to crossing orbits, they will 
physically collide with each other after a short time, as 
determined by the area of the planet subtended on the 
sky, 


47ra 2 m / a \ 2 {R v \ 

i “"~^M Tl ~ 10jTS (( mau) UJ 




V12days/ 

( 16 ) 

This expression assumes that the planets have acquired 
significant eccentricities and inclinations before impact, 
and that there is little gravitational focussing. Both as¬ 
sumptions are born out by our numerical simulations. 
In particular, the latter assumption is reasonble because 
at collision, their relative velocities can be higher than 
their surfa ce escape velocities (~ 10 — 20km/s, see, e.g. 


Wu & Lithwick 


2013), since the Keplerian orbital speed 


is of order 100 km/ s at the relevant distance. Such a 
collision will likely be destructive, as opposed to con- 


glome rative (Asphaug et al. 12006 Leinhardt & Stewart 


20121. Some planets may be lost arid massive debris may 


be produced. The former will increase the spacing within 
the system, while the latter will tend to cool down the or¬ 
bit and prevent further instabilities from happening. At 
the moment, we are incapable of predicting the outcome 
of such a collision, but speculate that this may explain 
the abundant low-multiple systems one observes in the 
Kepler field. 


4.2. Formation Environment 

We speculate briefly on the physical environment that 
can give rise to the compact spacing of Kepler systems. 
If our sculpting hypothesis proves correct, this is not re¬ 
stricted to only the high multiple systems, but likely ap¬ 
plies to all close-in planetary systems. 

The orbital spacings between terrestrial planets are 
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wide, with iVT ^ 30 — 40{^] These values have been repro- 


conglomerate out of a ses 

of smaller bodies (e.g., 

Cham- 

bers|2001; Kokubo et al. 

2006; O’Brien et al. 2006; Mor- 

ishima et al. 

2010). In tl 

lese simulations where there is 


the proto-planets can stir each other up to of order their 
surface escape velocities. Their orbits, being dynamically 
hot, require larger spacing to avoid further disru ption. 
This reasonin g is supported by simulations of Kominami 


& Ida (2002), which we discuss below. 

Starting from a population of Mars-sized proto¬ 
planets, but exerting an artificial dissipation on the 
plan et eccentricity with timescale Td amp , jKominami fe 


Ida (2002) showed that, terrestrial-like systems (A ~ 30J 
can be produced if Tdamp is relatively weak (more than 
a million times the orbital period), and closely-packed 
systems result if damping is strong. When Td amp ~ 10 5 
orbital time, in particular, the simulations produce sys¬ 
tems that are spaced between 8 to 13I?h - much like the 
initial population that we inferred in Fig. [H] At 0.1 AU, 
this corresponds to a damping time of Td amp ~ 3000 yrs. 
In Kepler systems, the primordial disk is likely more mas¬ 
sive than the terrestrial system studied by [Kominami & 


Ida (2002), so the requisite e-damping time may be even 


shorter. To put this timescale into context, the expected 
e-damping time from Type-I migration in a minimum 
mass solar nebula is Td amp ~ 3 yrs (Tanaka & Ward 


2004). 


So if our sculpting hypothesis is correct, we infer that 
the formation environment of the Kepler planets have 
to be moderately dissipative, with an eccentricity damp¬ 
ing timescale that is longer than the Type-I damping 
timescale, but much shorter than the million-year life¬ 
time of the nascent disks. 


4.3. Constraints on Orbital Shapes 

We can look at our results in a different way. The 
observed spacing could be used to constrain the values 
of eccentricity and mutual inclination in (and onl y in ) 
these high multiple systems. Inverting equation (14), 
adopting A' 5O (0,0) = 10.2 (required at 1 Gyrs) and 
Ji 5 o(cr e , cr irlc ) = 12.3 (observed at 1 Gyrs), we find the 
following independent limits, 


a e < 0.02 , and, 


< 0.084 radian 


(17) 


Although this constraint is obtained for our choice of 
mass-radius relationship, it likely applies even when the 
latter deviates from our choice. Similar constraints on 


orbital shapes were also obtained by (Mahajan & Wu 
20141 for a specific planetary system, Kepler-11. If cr. 


oy/ 2, as characteristic of 2-body relaxation in disks (Ida 
ct al. 1903). we would obtain a e < 0.02 and <j inc < 0.01 


0.0 U . These are indeed very cold dynamical systems. 

How do these constraints compare to what we have 
found via other channels? Most previous studies have 
not isolated high-multiple systems from the general cat¬ 
alogue and their results are therefore not directly com- 


inclinations have indicated low values (e.g. Fang & Mar- 

got 2012 

Tremaine & Dong 2012 

Fabrycky et al. 2014p. 


2 However, even this wide spacing, given the presence of giant 
planets and orbital eccentricities, still puts the Solar system at the 
edge of instability. 


There are only a few studies on planet eccentricities. 
Studying the distri bution of transit duratio ns in Kepler 


candidates have led Moorhead et al. (2011) to conclude 
that <r e ~ 0.2. However, this result is heavily plagued by 
uncertainties in stellar radius determinati ons an d is likely 


an over-estimate. TTY studies (Wu & Lithwick 2013 r 
Hadden & Lithwick 2014), on the other hand, have led to 


a much smaller eccentricity dispersion ( a e ~ 0.02), for a 
sub-population of planet pairs that are near mean-motion 
resonances. In fact, this small but non-zero dispersion 
partially motivates us in claiming that the observed sys¬ 
tems lie on the edge of instability. TTV studies also 
reported the existence of another sub-population where 
the eccentricities are substantially higher (tx e ~ 0.1), but 
these likely occur exclusively in low-multiple systems. 


5. CONCLUSION 

The architecture of Kepler multiple systems contains 
important pieces of clue for planet formation. In this 
work, we show that high-multiple Kepler systems have 
spacings that narrowly afford them dynamical stability 
in the giga-year timescale. This leads us to suggest that 
many more planetary systems were initially formed with 
even smaller spacing, but they have been since whittled 
down by dynamical instability, possibly producing some 
of the low-multiple Kepler systems we see today. If this 
is true, the final formation stage of these planets pro¬ 
ceed in highly dissipative environment, with eccentricity 
damping timescales much shorter than the disk lifetimes. 

Along the way, our numerical explorations have uncov¬ 
ered a number of interesting features in the stability of 
high-N systems. Consistent with previous investigations, 
we find that system stability time, quantified here as the 
time by which 50% of the initial systems remain stable, 
rises sharply with the spacing. When the spacings in a 
given system are not uniform but exhibit a spread, one 
can easily quantify the effect by adopting an effective 
spacing, K e g = A" mean — 0.5<r^. We give fitting formula 
for stability time when the planetary orbits are slightly 
eccentric and/or mutually inclined. Eccentricity, even at 
the level 1%, affects critical spacing significantly. As a 
result, we are able to place a stringent limit on the eccen¬ 
tricity of high-multiple Kepler systems, a e < 0.02. We 
find that systems that contain mean-motion-resonance 
pairs are preferentially destabilized. 

A large number of both theoretical and observational 
questions remain unanswered. What provides the strong 
eccentricity dissipation in planet formation? What is the 
origin of dynamical instability in high-multiple systems? 
What are the final outcomes in these dynamically unsta¬ 
ble systems? How do these outcomes compare against 
the observed systems? 
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